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Liquid  water  transport  in  the  diffusion  porous  layers  of  polymer  electrolyte  membrane  fuel  cells  (PEMFC) 
is  analyzed  as  a  process  of  quasi-static  invasion  from  multiple  interfacial  injection  sources.  From  pore 
network  simulations  based  on  a  new  version  of  the  invasion  percolation  algorithm  it  is  shown  that  a 
porous  layer  acts  as  a  two-phase  filter:  the  number  of  breakthrough  points  is  significantly  lower  that  the 
number  of  injection  points  owing  to  the  merging  of  liquid  paths  within  the  porous  layer.  The  number  of 
breakthrough  points  at  the  gas  diffusion  layer/gas  channel  interface  obtained  with  this  model  is  consistent 
with  the  available  experimental  observations. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Polymer  electrolyte  membrane  fuel  cells  (PEMFC)  are  regarded 
as  one  of  the  most  promising  devices  to  solve  the  energy  and  envi¬ 
ronmental  pollution  issues  [1].  As  discussed  in  many  works,  e.g. 
[1,2]  and  references  therein,  water  is  essential  for  the  operation 
of  PEMFC  and  the  so-called  water  management  problem  is  one  of 
the  key  issue  for  the  operation  of  PEMFC.  The  water  management 
problem  refers  to  the  fact  that  the  membrane  should  be  sufficiently 
hydrated  whereas  excessive  condensed  water  should  be  avoided 
in  other  parts  of  the  device,  such  as  the  catalyst  layer  (CL)  or  the 
gas  diffusion  layer  (GDL),  so  as  to  maintain  a  good  access  of  reac¬ 
tant  gas  to  reaction  sites.  A  popular  approach  to  this  problem  is 
to  rely  on  numerical  simulations  so  as  to  predict  and  to  analyse 
the  effects  of  various  parameters  on  pore  blockage  by  the  water. 
Most  of  these  studies  are  based  on  the  traditional  continuum  mod¬ 
els  describing  two-phase  flow  in  porous  media.  However,  there  are 
evidences  that  this  type  of  model  can  be  seriously  in  error,  notably 
owing  to  the  thin  nature  of  the  porous  layers  (lack  of  length  scale 
separation)  or  because  the  water  flow  regime  is  dominated  by  cap¬ 
illary  effects,  a  regime  which  cannot  be  described  using  continuum 
models,  see  [3]  and  references  therein  for  more  details.  As  consid¬ 
ered  in  several  previous  works  [3-6],  an  alternative  for  a  better 
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understanding  of  water  invasion  in  PEMFC  porous  layers  is  to  rely 
on  pore  network  models  and  associated  theories.  Since  water  inva¬ 
sion  is  dominated  by  capillary  effects  and  the  porous  layers  are 
generally  assumed  as  hydrophobic  (this  aspect  is  commented  at 
the  end  of  the  paper)  water  invasion  in  the  GDL  has  been  simulated 
using  the  invasion  percolation  (IP)  algorithm  [7]  since  it  is  well 
established  that  the  capillarity  dominated  displacement  of  a  wet¬ 
ting  fluid  (the  gas  phase  in  the  PEMFC  problem)  by  a  non-wetting 
one  (water  in  the  PEFMC  problem)  can  be  adequately  simulated 
with  the  IP  algorithm.  However,  as  pointed  out  in  [2]  and  [3],  there 
is  a  fundamental  difference  between  the  PEMFC  problem  and  the 
traditional  situations  involving  a  drainage  (drainage  refers  to  the 
process  of  displacing  a  wetting  fluid  by  a  non-wetting  one)  process 
at  very  low  capillary  numbers  in  porous  media  (the  capillary  num¬ 
ber  characterizes  the  competition  between  capillary  and  viscous 
forces,  the  later  can  be  neglected  for  a  sufficiently  small  capillary 
number).  In  the  traditional  use  of  the  IP  algorithm,  all  pores  at  the 
porous  medium  inlet  are  supposed  to  be  connected  to  a  reservoir 
and  the  invasion  is  simulated  increasing  steps  by  steps  very  gently 
the  pressure  of  the  non-wetting  fluid  in  the  reservoir.  In  the  PEMFC 
problem,  water  entering  a  porous  layer,  for  example  the  GDL,  comes 
from  an  adjacent  porous  layer,  the  micro-porous  layer  (MPL)  or  the 
CL,  and  not  from  a  reservoir  at  uniform  pressure.  Hence  the  injec¬ 
tion  conditions  in  the  PEMFC  problem  are  markedly  different  from 
the  traditional  injection  condition.  It  is  therefore  surmised  that  this 
has  a  great  impact  on  the  phase  distribution  and  therefore  the  pore 
blockage  by  the  water  [2,3].  A  well-known  fact  in  favour  of  this 
fundamental  difference  is  that  emergence  of  several  liquid  water 
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droplets  is  generally  observed  at  the  GDL/gas  channel  interface 
[8-10].  This  cannot  be  predicted  when  the  IP  algorithm  is  used  with 
the  traditional  “connected  reservoir”  boundary  condition  since  only 
one  breakthrough  droplet  is  predicted  in  this  case.  If  condensation 
effects  and  viscous  effects  can  be  neglected,  one  obvious  possibility 
for  predicting  multiple  breakthrough  points  with  the  IP  algorithm 
is  to  consider  that  water  enters  into  a  porous  layer  through  mul¬ 
tiple  independent  inlet  injection  points  [2,11].  One  argument  in 
favour  of  this  scenario  is  that  the  MPL  pores  are  approximately 
about  100  times  smaller  than  those  in  the  GDL  and  that  the  size 
of  pores  between  the  agglomerates  in  the  CL  is  in  turn  at  least  10 
times  smaller  than  the  average  pore  size  in  the  MPL.  Because  of  this 
hierarchy  in  pore  sizes  liquid  water  produced  in  the  CL  can  enter  the 
MPL  from  many  pores.  This  leads  to  the  possibility  of  formation  of 
several  liquid  paths  into  the  MPL.  If  there  is  no  complete  merging  of 
those  liquid  paths  inside  the  MPL,  this  leads  in  turn  to  the  possibil¬ 
ity  of  several  injection  points  into  the  GDL  and  eventually,  if  again 
there  is  no  complete  merging  of  liquid  paths  within  the  GDL,  to 
the  possible  emergence  of  several  droplets  at  the  GDL/gas  channel 
interface.  In  this  article,  we  explore  the  relevance  of  this  scenario 
from  IP  simulations  with  a  new  inlet  boundary  condition  consid¬ 
ering  that  liquid  water  can  enter  each  porous  layer  from  multiple 
injection  points. 


2.  IP  algorithm  with  multiple  inlet  injections 

As  in  a  previous  work  [3],  we  consider  that  the  size  of  a  unit  cell 
of  the  system  is  fixed  by  the  mean  distance  L  between  two  chan¬ 
nels  in  the  bipolar  plate  (assuming  for  simplicity  a  system  of  parallel 
channels)  with  L  ^  2  mm.  This  leads  to  study  the  liquid  invasion  in 
porous  systems  of  size  Lx  Lx  l  where  L  is  the  porous  layer  thick¬ 
ness.  As  discussed  in  [3],  a  representative  mean  distance  a  between 
two  pores  in  the  GDL  can  be  estimated  as  a  &  50  |xm  whereas  l  is 
typically  in  the  range  [  1 70-400]  |xm.  If  for  simplicity  we  do  not  take 
into  account  the  GDL  microstructure  anisotropy  (this  aspect  will 
be  considered  in  a  future  work),  this  leads  to  the  consideration  of 
40  x  40  x  N  pore  networks  with  N  varying  in  the  range  [4, 40]  (with 
a  lattice  spacing  a  of  50  |jim,  N= 4  corresponds  for  example  to  a  GDL 
thickness  of  200  |xm).  The  MPL  is  made  of  particles  and  has  pore 
sizes  of  the  order  of  0.5  |xm  whereas  its  thickness  is  typically  in  the 
range  10-20  pan.  As  a  reasonable  guess  of  mean  distance  between 
two  pores  in  MPL,  we  take  2  [xm.  This  leads  to  1 000  x  1 000  x  M  pore 
networks  with  M  varying  in  the  range  [5,10]  (with  a  lattice  spac¬ 
ing  a  of  2  pm,  M  =  5  corresponds  for  example  to  a  MPL  thickness  of 
10  pm).  As  in  [3],  the  simulations  are  performed  over  simple  cubic 
networks  of  pores  connected  by  throats,  with  the  throat  and  pore 
sizes  randomly  distributed  over  the  network  as  explained  in  [3]. 

Since  water  is  produced  in  the  agglomerates  of  the  CL,  i.e.  poten¬ 
tially  everywhere  in  the  CL,  and  since  the  size  of  secondary  pores 
in  the  CL  (the  secondary  pores  are  the  pores  between  the  agglom¬ 
erates)  is  about  10  times  smaller  than  the  pore  size  in  the  MPL, 
we  assume  here  for  simplicity  that  water  enters  the  MPL  in  liq¬ 
uid  form  at  the  CL/MPL  interface  through  all  pores  of  MPL  located 
at  the  CL/MPL  interface.  Hence  we  assume  that  the  number  nmi  of 
pores  at  the  CL/MPL  interface  where  liquid  water  penetrates  into 
the  MPL  is  nmi  =  1000  x  1000.  In  fact  we  expect  that  nmi  varies  with 
the  current  density  but  this  would  be  explored  in  a  future  work. 
To  determine  the  number  of  breakthrough  points  nmb  at  the  MPL 
outlet  (MPL/GDL  interface),  we  use  the  following  IP  algorithm.  We 
determine  the  flow  path  using  the  standard  IP  algorithm  without 
trapping  [7]  starting  with  only  one  injection  point  (the  first  pore  in 
the  first  pores  row  at  the  MPL  inlet).  Then  we  repeat  the  simula¬ 
tion  starting  from  the  second  pore  in  the  first  row.  We  stops  this 
second  simulation  either  when  the  flow  path  generated  from  this 
second  injection  point  merges  into  the  flow  path  associated  with 


Fig.  1-  Evolution  of  average  penetration  (J)  as  a  function  of  layer  thickness  M  (MPL) 
or  N  (GDL)  measured  in  lattice  spacing  unit.  The  results  for  the  MPL  1 000  x  1 000  x  M 
networks  (black  crosses)  and  for  the  GDL  40  x  40  x  N  networks  (red  circles)  collapse 
in  a  single  curve  except  for  the  larger  values  of  N  for  which  there  are  less  than  10 
breakthrough  points.  The  inset  shows  a  log-log  plot  of  the  same  data  together  with 
a  best  linear  fit  showing  that  (j)  ~N-1-87  for  N  or  M  >  6. 

the  first  inlet  injection  pore  or  at  breakthrough,  i.e.  when  the  liquid 
injected  from  the  second  inlet  pore  reaches  the  MPL/GDL  interface 
through  a  path  independent  from  the  path  connected  to  the  first 
injection  point.  Then  we  repeat  this  procedure  starting  successively 
from  all  the  other  pores  at  the  CL/MPL  interface.  Spatial  periodicity 
boundary  conditions  are  imposed  on  the  lateral  sides  of  network. 

3.  Results 

So  as  to  take  into  account  the  statistical  fluctuations  inherent  to 
a  IP  process  on  a  finite  size  structure,  the  simulations  are  performed 
over  200  realizations,  unless  otherwise  mentioned,  of  the  cubic  lat¬ 
tice  for  each  size  M  considered.  The  results  are  analyzed  in  terms  of 
the  “penetration” /(by  analogy  with  particles  filtration  problems), 
which  is  defined  as  f=nmblnmi.  Thus  /  represents  the  fraction  of 
pores  at  the  porous  medium  outlet  where  a  breakthrough  occurs, 
or  alternatively  the  probability  that  a  pore  at  the  porous  medium 
outlet  is  a  breakthrough  pore.  The  evolution  of  (J),  where  (j)  is  the 
penetration  averaged  over  the  number  of  realizations  considered, 
is  shown  in  Fig.  1  as  a  function  of  M.  As  can  be  seen  from  Fig.  1, 
(j)  drops  rapidly  with  M.  The  MPL  acts  in  fact  as  a  quite  efficient 
two-phase  filter  since  the  fraction  if)  of  breakthrough  pores  at  the 
MPL  outlet  is  significantly  lower  than  1 ,  even  for  M  as  small  as  1  for 
which  (f)  ^  0.3.  As  shown  in  the  inset  in  Fig.  1,  the  variation  of  (j) 
with  M  is  consistent  with  a  power  law  of  the  form  (j)  =p(£/a)-1,87 
with p  ^  0.882,  except  for  M  =  L\a<  6.  It  is  expected  that  this  scaling 
holds  only  over  a  limited  range  of  M  for  M> 6  since  the  number  of 
breakthrough  pore  nmb  cannot  be  lower  that  1 .  This  is  illustrated 
further  below  when  the  GDL  smaller  network  is  considered.  This 
indicates  that  the  merging  of  liquid  flow  paths  eventually  leads  to 
leave  many  pores  free  of  water  and  therefore  accessible  to  gas  even 
for  the  extreme  case  considered  here  where  all  pores  at  the  inlet  are 
injection  points.  Owing  to  the  randomness  in  the  throats  size,  the 
number  of  breakthrough  pores  nmb ,  and  therefore  the  penetration 
/,  are  random  variables.  As  will  be  illustrated  below,  the  probability 
density  function  (p.d.f.)  of  the  number  of  breakthrough  points  is 
Gaussian.  The  values  of  ( nmb )  and  of  the  standard  deviation  on  of 
nmb  obtained  from  our  simulation  are  summarized  in  Table  1. 
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Table  1 

Evolution  of  average  number  of  breakthrough  points  {nmb)  at  MPL  outlet  and  stan¬ 
dard  deviation  of  nmh  as  a  function  of  MPL  thickness  (measured  in  lattice  spacing). 


M 

(nmb) 

On 

1 

297320.0 

443.0 

2 

142789.0 

290.5 

3 

83945.8 

200.2 

4 

55020.0 

166.5 

5 

38744.2 

127.8 

6 

28705.4 

108.8 

7 

22096.5 

110.7 

8 

17504.8 

87.4 

9 

14200.7 

87.2 

10 

11741.1 

72.91 

11 

9877.0 

63.25 

12 

8404.0 

56.97 

13 

7249.18 

57.71 

14 

6307.5 

50.02 

15 

5542.54 

49.16 

16 

4907.28 

52.64 

17 

4380.52 

41.33 

18 

3927.34 

42.31 

19 

3538.76 

38.94 

20 

3211.23 

35.24 

Table  2 

Evolution  of  average  number  of  breakthrough  points  (ngb)  at  GDL  outlet,  standard 
deviation  of  ngb  and  breakthrough  number  density  as  a  function  of  GDL  thickness 
(measured  in  lattice  spacing). 


N 

(ftgb) 

Cfn 

Breakthrough  number 
densities  (pores  per  mm2) 

1 

475.4 

17.43 

119 

2 

228.4 

9.91 

57 

3 

134.8 

8.158 

33 

4 

88.07 

6.447 

22 

5 

61.67 

5.604 

15 

6 

45.87 

4.645 

11 

7 

35.76 

4.026 

8 

8 

28.12 

3.417 

7 

9 

22.93 

3.333 

5 

10 

18.67 

2.98 

4 

11 

15.82 

2.794 

3 

12 

13.63 

2.521 

3 

13 

11.61 

2.441 

3 

14 

10.23 

2.086 

2 

15 

8.72 

2.018 

2 

16 

7.88 

1.961 

1.75 

17 

7.1 

1.944 

1.75 

18 

6.415 

1.686 

1.5 

19 

5.805 

1.666 

1.25 

20 

5.285 

1.563 

1.25 

21 

4.725 

1.562 

1 

22 

4.075 

1.375 

1. 

23 

3.84 

1.365 

0.75 

24 

3.665 

1.274 

0.75 

25 

3.265 

1.29 

0.75 

26 

3.08 

1.214 

0.75 

27 

2.87 

1.124 

0.50 

28 

2.76 

1.115 

0.50 

29 

2.665 

1.001 

0.50 

30 

2.42 

1.115 

0.50 

31 

2.32 

0.9042 

0.50 

32 

2.175 

0.8913 

0.50 

33 

2 

0.8544 

0.50 

34 

1.855 

0.7576 

0.25 

35 

1.78 

0.7427 

0.25 

36 

1.795 

0.7956 

0.25 

37 

1.645 

0.7609 

0.25 

38 

1.645 

0.7064 

0.25 

39 

1.55 

0.6225 

0.25 

40 

1.395 

0.5648 

0.25 

Fig.  2.  Semi-log  plot  of  standard  deviation  a  of  penetration/as  a  function  of  average 
penetration  if)  (the  red  solid  line  with  crosses  corresponds  to  the  GDL  40  x  40  x  N 
networks  whereas  the  black  solid  line  with  circles  is  for  the  MPL  1000  x  1000  x  M 
networks).  The  inset  shows  the  p.d.f.  (probability  density  function)  of  ngb  for  N=4 
obtained  from  5000  realizations  of  the  GDL  network. 


To  determine  the  number  of  breakthrough  points  ngb  at  the  GDL 
outlet,  we  use  the  same  IP  algorithm  as  before.  A  possible  essential 
difference  is  that  not  all  pores  at  the  GDL  inlet  (MPL/GDL  interface) 
are  necessarily  injection  points  since  water  exits  the  MPL  through 
a  limited  number  of  breakthrough  points,  see  Table  1.  As  noted 
before,  the  mean  distance  between  two  pores  is  of  the  order  of  2  fxrn 
in  the  MPL  and  50  p,m  in  the  GDL.  The  simulations  show  that  the 
breakthrough  points  are  randomly  distributed  over  the  MPL  outlet 
surface.  Hence  when  the  number  of  breakthrough  points  nmb  at  the 
MPL  outlet  is  lower  than  40  x  40  =  1600,  a  reasonable  assumption 
is  that  there  is  nmb  injection  pores  into  the  GDL  since  it  is  likely  that 
each  breakthrough  points  from  the  MPL  is  connected  to  a  different 
pore  of  the  GDL  at  the  MPL/GDL  interface,  that  is  =  nmb.  When 
nmb  is  greater  than  1600,  then  several  MPL  breakthrough  pores  can 
be  connected  to  the  same  interfacial  pore  of  the  GDL  and  therefore 
rigi  =  1 600  since  ngi  cannot  be  greater  than  1 600.  The  values  reported 
in  Table  1  show  that  nmb  is  greater  than  1600  for  the  range  of  MPL 
thickness  considered  and  when  all  MPL  inlet  pores  are  injection 
pores.  We  therefore  explore  first  the  case  ngi  =  1600. 

The  evolution  of  average  penetration  (f)  for  the  GDL  network  is 
shown  in  Fig.  1  (red  line  with  circles)  while  the  inset  in  Fig.  1  shows 
the  same  data  in  log-log  coordinates.  As  can  be  seen  there  is  a  per¬ 
fect  collapse  of  the  data  obtained  from  the  MPL  1000  x  1000  x  M 
network  and  from  the  GDL  40  x  40  x  N  network  except  for  the 
larger  values  of  N  for  which  the  number  of  breakthrough  points 
becomes  lower  than  10,  that  is  too  small  for  observing  a  quasi- 
continuous  evolution  of  if)  with  N.  Hence  the  staircase  evolution  of 
(J)  reflects  the  fact  that  there  are  only  a  few  breakthrough  points 
for  the  larger  N  (see  Table  2).  As  mentioned  before  and  shown  in 
the  inset  in  Fig.  2,  the  p.d.f  of  the  number  of  breakthrough  pores  is 
Gaussian.  Contrary  to  the  evolution  of  if),  which  is  independent  of 
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Fig.  3.  Water  saturation  profiles  along  the  GDL  for  N=  10  for  various  average  num¬ 
ber  of  inlet  injection  points  ng,  (1600,  800,  400,  200, 100,  50)  randomly  distributed 
over  the  GDL  inlet  face.  Each  profile  is  an  average  over  200  realizations  of  the  GDL 
network.  The  saturation  fluctuations  are  small,  so  that  a  variation  of  one  standard 
deviation  about  the  average  values  would  not  be  discernible  in  the  figure. 

L/a,  there  is  a  significant  effect  of  network  size  L\a  on  the  standard 
deviation  off  This  is  shown  in  Fig.  2.  As  can  be  seen,  ctflf)  decreases 
with  network  size  with  ay/  (J)  in  the  range  [  1 0-3 , 1 0-2  ]  for  L\a  =  1 000 
and  (J f  I  (J)  two  orders  of  magnitude  greater  for  L/a  =  40.  This  aspect 
will  be  explored  in  more  detail  in  a  forthcoming  work. 

The  breakthrough  density  dn  (number  of  breakthroughs  per  unit 
area)  at  the  GDL  outlet  (GDL/channel  interface)  obtained  in  our  sim¬ 
ulations  is  reported  in  Table  2.  As  can  be  seen,  (dn)  is  in  the  range 
[4-22]  pores  per  mm2  for  porous  layer  thicknesses  representative 
of  those  of  GDLs.  This  is  in  very  good  agreement  with  the  observa¬ 
tions  of  droplet  formations  at  the  GDL/channel  interface  reported 
in  [10]. 

As  a  final  result,  Fig.  3  shows  the  liquid  saturation  profiles  along 
the  GDL  thickness  for  N=  10.  These  profiles  are  obtained  by  deter¬ 
mining  the  saturation  in  successive  slices  of  the  network  along  the 
GDL  thickness.  As  stated  before,  it  can  be  surmised  that  the  num¬ 
ber  of  injection  points  varies  with  the  current  density.  So  as  to  have 
a  first  insight  into  the  impact  of  the  number  of  injection  points, 
the  profiles  corresponding  to  various  average  numbers  of  injection 
points  are  shown  in  Fig.  3.  Interestingly,  the  profile  shape  changes 
from  concave  to  convex  as  the  number  of  injection  points  decreases. 
As  discussed  for  example  in  [3,1 1  ],  the  shape  of  the  saturation  pro¬ 
file  is  an  important  indication  as  regards  the  analysis  of  two-phase 
flow  in  GDL.  Lienee  our  results  can  explain  why  profiles  deduced 
from  in  situ  experiments  [12,13]  are  not  concave.  Also,  it  can  be  seen 
from  Fig.  3  that  the  saturations  along  the  GDL  thickness  are  here 
significantly  greater  than  when  the  traditional  “connected  reser¬ 
voir”  boundary  condition  is  used  (see  for  example  Fig.  14  in  [3]), 
especially  in  the  second  half  of  profiles.Hence  the  traditional  use 
of  the  IP  model  with  the  “connected  reservoir”  boundary  condition 


leads  to  seriously  underestimate  the  pore  blockage  by  the  water. 
The  inlet  independent  injection  points  boundary  condition  appears 
as  much  more  appropriate  for  analyzing  two-phase  flows  in  a  GDL 
as  well  as  in  the  MPL. 

4.  Conclusion 

We  have  explored  a  possible  scenario  of  droplet  generations  at 
the  GDL/gas  channel  interface  assuming  that  water  invasion  into 
the  GDL  and  the  MPL  could  be  described  using  the  IP  algorithm.  Con¬ 
trary  to  classical  IP  simulations,  the  present  simulations  are  based 
on  the  concept  that  liquid  water  enters  the  porous  layer  through 
multiple  independent  injection  sources  at  the  inlet.  The  results  are 
consistent  with  the  available  data  on  the  density  of  droplets  form¬ 
ing  at  the  GDL/gas  channel  interface  as  well  as  with  the  saturation 
profile  shape  observed  in  in  situ  experiments.  The  algorithm  pro¬ 
poses  in  this  article  provides  one  tool  to  better  analyze  the  water 
management  problem  in  PEMFC.  Flowever,  the  droplets  generation 
process  at  the  surface  of  a  porous  layer  described  in  this  article 
is  probably  not  the  only  one  taking  place  in  PEMFC.  Condensa¬ 
tion  processes  (more  likely  however  on  the  anode  side  than  on 
the  cathode  side  [2])  and  viscous  effects  [11]  (neglected  in  the 
present  effort)  could  also  play  a  role  in  the  droplet  formations. 
Also,  we  have  assumed  perfectly  hydrophobic  porous  layers  and 
isotropic  microstructures.  As  considered  for  example  in  [14],  there 
are  evidences  that  the  GDL  is  only  partially  hydrophobic  and  it 
is  well  known  that  the  GDL  microstructure  is  not  isotropic.  The 
effects  of  a  mixed  wettability  and  anisotropy  on  droplets  genera¬ 
tions  would  deserve  to  be  explored  using  an  approach  similar  to 
the  one  described  in  this  article. 
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